Chapter 4 Diversity analysis
genome_counts_filt <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
select(where(~ sum(.) > 0)) %>%
select(-AC85) %>%
rownames_to_column(., "genome")
genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt$genome)
table <- genome_counts_filt %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character,as.numeric) %>%
rownames_to_column(., "sample")
master_table <- sample_metadata %>%
mutate(sample=Tube_code) %>%
mutate(Tube_code=str_remove_all(Tube_code, "_a")) %>%
filter(Tube_code %in% table$sample) %>%
mutate(treatment = sub("^\\d+_", "", time_point))%>%
left_join(., table, by=join_by("Tube_code" =="sample"))
sample_metadata <- master_table %>%
select(1:13)
genome_counts_filt <- master_table %>%
select(12,14:323) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character,as.numeric) %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
dplyr::select(where(~ !all(. == 0)))%>%
rownames_to_column(., "genome")
genome_counts_filtering <- master_table %>%
select(12,14:323) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character,as.numeric) %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0)%>%
dplyr::select(where(~ !all(. == 0)))
genome_tree <- keep.tip(genome_tree, tip=genome_counts_filt$genome)
#match_data(genome_counts_filtering,genome_tree)
genome_metadata <- genome_metadata %>%
filter(genome %in% genome_counts_filt$genome)4.1 Calculate diversities
4.1.1 Alpha diversity
# Calculate Hill numbers
richness <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 0) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(richness = 1) %>%
rownames_to_column(var = "sample")
neutral <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(neutral = 1) %>%
rownames_to_column(var = "sample")
phylogenetic <- genome_counts_filt %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, tree = genome_tree) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(phylogenetic = 1) %>%
rownames_to_column(var = "sample")
genome_gifts1 <- genome_gifts[genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]
genome_counts_filt1 <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts1),]
genome_counts_filt1 <- genome_counts_filt1 %>%
remove_rownames() %>%
column_to_rownames(var = "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
rownames_to_column(., "genome")
dist <- genome_gifts1 %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
functional <- genome_counts_filt1 %>%
filter(genome %in% rownames(dist)) %>%
column_to_rownames(var = "genome") %>%
dplyr::select(where(~ !all(. == 0))) %>%
hilldiv(., q = 1, dist = dist) %>%
t() %>%
as.data.frame() %>%
dplyr::rename(functional = 1) %>%
rownames_to_column(var = "sample") %>%
mutate(functional = if_else(is.nan(functional), 1, functional))
# Merge all metrics
alpha_div <- richness %>%
full_join(neutral, by = join_by(sample == sample)) %>%
full_join(phylogenetic, by = join_by(sample == sample))%>%
full_join(functional, by = join_by(sample == sample))4.1.2 Beta diversity
beta_q0n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 0)
beta_q1n <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1)
beta_q1p <- genome_counts_filt %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, tree = genome_tree)
genome_gifts1 <- genome_gifts[genome_counts_filt$genome,]
genome_gifts1 <- genome_gifts1[, colSums(genome_gifts1 != 0) > 0]
genome_counts_filt1 <- genome_counts_filt[genome_counts_filt$genome %in% rownames(genome_gifts1),]
dist <- genome_gifts1 %>%
to.elements(., GIFT_db) %>%
traits2dist(., method = "gower")
beta_q1f <- genome_counts_filt1 %>%
column_to_rownames(., "genome") %>%
filter(rowSums(. != 0, na.rm = TRUE) > 0) %>%
select_if(~!all(. == 0)) %>%
hillpair(., q = 1, dist = dist)4.2 Does captivity change the microbial community?
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "0_Wild") ) 4.2.1 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "0_Wild") ) %>%
ggplot(aes(y = value, x = time_point, color=time_point, fill=time_point)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(values=treatment_colors)+
scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")Richness
Modelq0GLMMNB <- glmer.nb(richness ~ time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(7.7073) ( log )
Formula: richness ~ time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
341.1 347.3 -166.5 333.1 31
Scaled residuals:
Min 1Q Median 3Q Max
-1.88922 -0.29614 0.04658 0.45758 1.21929
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0.1706 0.413
Number of obs: 35, groups: individual, 18
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 4.1983 0.1330 31.564 < 2e-16 ***
time_point1_Acclimation -0.4405 0.1344 -3.279 0.00104 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
tm_pnt1_Acc -0.456
Neutral
Model_neutral <- lme(fixed = neutral ~ time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
258.7499 264.7359 -125.3749
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 8.169818 7.178627
Fixed effects: neutral ~ time_point
Value Std.Error DF t-value p-value
(Intercept) 36.65610 2.563403 17 14.299780 0
time_point1_Acclimation -16.81743 2.447304 16 -6.871818 0
Correlation:
(Intr)
time_point1_Acclimation -0.456
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.52249470 -0.58253897 -0.03654165 0.58990386 1.40500351
Number of Observations: 35
Number of Groups: 18
Phylogenetic
Model_phylo <- lme(fixed = phylogenetic ~ time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
129.7532 135.7392 -60.8766
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.6806603 1.248658
Fixed effects: phylogenetic ~ time_point
Value Std.Error DF t-value p-value
(Intercept) 5.221589 0.3351985 17 15.577603 0.0000
time_point1_Acclimation 0.135310 0.4236755 16 0.319371 0.7536
Correlation:
(Intr)
time_point1_Acclimation -0.61
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.18114503 -0.47484996 0.08008291 0.37909497 1.82560478
Number of Observations: 35
Number of Groups: 18
Functional
Model_funct <- lme(fixed = functional ~ time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_funct)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
-52.38562 -46.39959 30.19281
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.03084103 0.08367718
Fixed effects: functional ~ time_point
Value Std.Error DF t-value p-value
(Intercept) 1.5154110 0.02101988 17 72.09416 0.0000
time_point1_Acclimation -0.0307069 0.02834791 16 -1.08322 0.2948
Correlation:
(Intr)
time_point1_Acclimation -0.653
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.20739889 -0.40096259 0.05458715 0.62964551 1.42654682
Number of Observations: 35
Number of Groups: 18
4.2.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "0_Wild")) %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "0_Wild"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.07277 0.072771 7.2819 999 0.012 *
Residuals 33 0.32979 0.009993
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(richness ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.5295984 | 0.05492497 | 1.917862 | 0.001 |
| Residual | 33 | 9.1126169 | 0.94507503 | NA | NA |
| Total | 34 | 9.6422153 | 1.00000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.02764 0.027637 2.2078 999 0.128
Residuals 33 0.41309 0.012518
adonis2(neutral ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.7470135 | 0.08504426 | 3.067318 | 0.001 |
| Residual | 33 | 8.0368067 | 0.91495574 | NA | NA |
| Total | 34 | 8.7838201 | 1.00000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04119 0.041187 2.897 999 0.089 .
Residuals 33 0.46917 0.014217
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(phylo ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.2075719 | 0.1133094 | 4.217042 | 0.001 |
| Residual | 33 | 1.6243310 | 0.8866906 | NA | NA |
| Total | 34 | 1.8319029 | 1.0000000 | NA | NA |
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00017 0.0001696 0.0078 999 0.931
Residuals 33 0.71884 0.0217831
adonis2(func ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(func))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(func))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.0009227428 | 0.0009988421 | 0.03299475 | 0.615 |
| Residual | 33 | 0.9228897312 | 0.9990011579 | NA | NA |
| Total | 34 | 0.9238124740 | 1.0000000000 | NA | NA |
4.2.3 Structural zeros
acclimation_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point == "1_Acclimation") %>%
dplyr::select(sample) %>%
pull()
wild_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point == "0_Wild") %>%
dplyr::select(sample) %>% pull()
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",subset_meta$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>%
mutate(all_zeros_wild= all(c_across(all_of(wild_samples)) == 0)) %>%
mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>%
mutate(average_wild = mean(c_across(all_of(wild_samples)), na.rm = TRUE)) %>%
filter(all_zeros_acclimation == TRUE || all_zeros_wild==TRUE) %>%
mutate(present = case_when(
all_zeros_acclimation & !all_zeros_wild ~ "wild",
!all_zeros_acclimation & all_zeros_wild ~ "acclimation",
!all_zeros_acclimation & !all_zeros_wild ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "acclimation", average_acclimation, average_wild)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Acclimation
structural_zeros %>%
filter(present=="acclimation") %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) # A tibble: 6 × 2
# Rowwise:
phylum sample_count
<chr> <int>
1 p__Pseudomonadota 5
2 p__Bacillota 3
3 p__Bacillota_A 2
4 p__Bacteroidota 2
5 p__Actinomycetota 1
6 p__Chlamydiota 1
Wild
# A tibble: 11 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_19:bin_000031 p__Bacillota_A c__Clostridia o__Lachnospirales f__Lachnospiraceae g__Lachnotalea s__
2 AH1_2nd_7:bin_000023 p__Bacillota_A c__Clostridia o__Lachnospirales f__ g__ s__
3 AH1_2nd_20:bin_000025 p__Bacillota_A c__Clostridia o__Lachnospirales f__CAG-274 g__ s__
4 AH1_2nd_6:bin_000024 p__Spirochaetota c__Brevinematia o__Brevinematales f__Brevinemataceae g__Brevinema s__
5 AH1_2nd_10:bin_000097 p__Bacteroidota c__Bacteroidia o__Bacteroidales f__Tannerellaceae g__Parabacteroides s__
6 LI1_2nd_7:bin_000020 p__Bacillota_B c__Peptococcia o__Peptococcales f__Peptococcaceae g__ s__
7 LI1_2nd_7:bin_000001 p__Bacillota_A c__Clostridia o__Clostridiales f__Clostridiaceae g__Sarcina s__
8 LI1_2nd_7:bin_000083 p__Bacillota_A c__Clostridia o__Oscillospirales f__Oscillospiraceae g__Pseudoflavonifractor s__
9 AH1_2nd_7:bin_000055 p__Bacillota c__Bacilli o__Mycoplasmatales f__Mycoplasmoidaceae g__Ureaplasma s__
10 LI1_2nd_1:bin_000052 p__Bacillota_A c__Clostridia o__TANB77 f__CAG-508 g__RGIG8482 s__
11 AH1_2nd_6:bin_000060 p__Bacillota_A c__Clostridia o__Oscillospirales f__Oscillospiraceae g__Intestinimonas s__
Community elements differences in structural zeros
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)# A tibble: 35 × 3
trait p_value p_adjust
<chr> <dbl> <dbl>
1 B0214 0.000261 0.00340
2 B0219 0.00164 0.0111
3 B0302 0.000198 0.00284
4 B0703 0.00306 0.0182
5 B0709 0.000103 0.00284
6 B0712 0.000551 0.00561
7 B0801 0.000198 0.00284
8 B1012 0.00136 0.0102
9 B1028 0.00000837 0.000599
10 D0102 0.00131 0.0102
# ℹ 25 more rows
4.2.4 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "0_Wild") )%>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "time_point",
rand_formula = "(1|individual)",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output$res %>%
dplyr::select(taxon, lfc_time_point1_Acclimation, p_time_point1_Acclimation) %>%
filter(p_time_point1_Acclimation < 0.05) %>%
dplyr::arrange(p_time_point1_Acclimation) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_time_point1_Acclimation)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_point1_Acclimation, y=forcats::fct_reorder(genome,lfc_time_point1_Acclimation), fill=phylum)) + #forcats::fct_rev()
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
# coord_flip()+
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))Phyla of the significant MAGs
ancombc_rand_table_mag%>%
filter(lfc_time_point1_Acclimation<0) %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) phylum sample_count
1 Bacillota_A 16
2 Bacteroidota 9
3 Bacillota_C 1
phylum genus species
1 Bacillota_A MGBC140009
2 Bacillota_A
3 Bacillota_A
4 Bacillota_A Negativibacillus
5 Bacteroidota Parabacteroides
6 Bacteroidota Bacteroides
7 Bacillota_A Hungatella
8 Bacteroidota Bacteroides
9 Bacteroidota Parabacteroides
10 Bacteroidota Parabacteroides
11 Bacteroidota Bacteroides
12 Bacillota_A
13 Bacteroidota Bacteroides
14 Bacillota_A
15 Bacillota_A Copromonas
16 Bacillota_A Enterocloster
17 Bacillota_A Velocimicrobium
18 Bacillota_A CAG-95
19 Bacillota_C
20 Bacteroidota Bacteroides
21 Bacillota_A
22 Bacillota_A Copromonas
23 Bacteroidota Bacteroides Bacteroides thetaiotaomicron
24 Bacillota_A Pseudoflavonifractor
25 Bacillota_A Enterocloster
26 Bacillota_A JALFVM01
ancombc_rand_table_mag%>%
filter(lfc_time_point1_Acclimation<0) %>%
count(genus, name = "sample_count") %>%
arrange(desc(sample_count)) genus sample_count
1 6
2 Bacteroides 6
3 Parabacteroides 3
4 Copromonas 2
5 Enterocloster 2
6 CAG-95 1
7 Hungatella 1
8 JALFVM01 1
9 MGBC140009 1
10 Negativibacillus 1
11 Pseudoflavonifractor 1
12 Velocimicrobium 1
4.2.5 Differences in the functional capacity
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Wild"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.2.5.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.331 0.0320
2 Wild 0.345 0.0234
MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94347, p-value = 0.07144
Wilcoxon rank sum exact test
data: value by treatment
W = 114, p-value = 0.2066
alternative hypothesis: true location shift is not equal to 0
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-Wild)%>%
mutate(group_color = ifelse(Difference <0, "Wild","Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=treatment_colors) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")4.2.5.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.329 0.0320
2 Wild 0.346 0.0195
MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94492, p-value = 0.07913
Wilcoxon rank sum exact test
data: value by treatment
W = 100, p-value = 0.08313
alternative hypothesis: true location shift is not equal to 0
function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))| Code_function | Acclimation | Wild | Function | Difference | treatment |
|---|---|---|---|---|---|
| B03 | 0.1097356 | 0.1256091 | Amino acid derivative biosynthesis | -0.0158735 | Wild |
| D03 | 0.2921075 | 0.3442827 | Sugar degradation | -0.0521752 | Wild |
4.3 Does the antibiotic administration change the microbial community?
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics") ) 4.3.1 Alpha diversity
alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics") ) %>%
ggplot(aes(y = value, x = time_point, color=time_point, fill=time_point)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(values=treatment_colors)+
scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")Richness
Modelq0GLMMNB <- glmer.nb(neutral ~ time_point+(1|individual), data = alpha_div_meta)
summary(Modelq0GLMMNB)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(4.0665) ( log )
Formula: neutral ~ time_point + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
223.0 228.8 -107.5 215.0 28
Scaled residuals:
Min 1Q Median 3Q Max
-1.4702 -0.6482 -0.1777 0.6405 2.6346
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 0.07573 0.2752
Number of obs: 32, groups: individual, 18
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.9069 0.1546 18.800 < 2e-16 ***
time_point2_Antibiotics -0.9393 0.2131 -4.407 1.05e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
tm_pnt2_Ant -0.538
Neutral
Model_neutral <- lme(fixed = neutral ~ time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
224.1303 229.7351 -108.0652
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 3.198228 7.479481
Fixed effects: neutral ~ time_point
Value Std.Error DF t-value p-value
(Intercept) 19.11841 1.971629 17 9.696759 0e+00
time_point2_Antibiotics -11.46314 2.675428 13 -4.284600 9e-04
Correlation:
(Intr)
time_point2_Antibiotics -0.63
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.86893024 -0.51778639 -0.08383536 0.57135437 1.94021414
Number of Observations: 32
Number of Groups: 18
Phylogenetic
Model_phylo <- lme(fixed = phylogenetic ~ time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
137.9324 143.5371 -64.96618
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 1.019275 1.66934
Fixed effects: phylogenetic ~ time_point
Value Std.Error DF t-value p-value
(Intercept) 5.400123 0.4734177 17 11.406677 0.000
time_point2_Antibiotics -0.791476 0.6015713 13 -1.315681 0.211
Correlation:
(Intr)
time_point2_Antibiotics -0.587
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.65991143 -0.50025273 -0.03029114 0.42477766 1.61395665
Number of Observations: 32
Number of Groups: 18
Functional
Model_funct <- lme(fixed = functional ~ time_point, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_funct)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
-15.05258 -9.447793 11.52629
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 1.585329e-06 0.1502428
Fixed effects: functional ~ time_point
Value Std.Error DF t-value p-value
(Intercept) 1.4848858 0.03643923 17 40.74965 0.000
time_point2_Antibiotics -0.0341744 0.05322290 13 -0.64210 0.532
Correlation:
(Intr)
time_point2_Antibiotics -0.685
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.2441309 -0.5590927 0.2912519 0.6141833 2.0349226
Number of Observations: 32
Number of Groups: 18
4.3.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics")) %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.013687 0.0136873 2.0494 999 0.162
Residuals 30 0.200356 0.0066785
adonis2(richness ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.580367 | 0.1313608 | 4.536779 | 0.001 |
| Residual | 30 | 10.450370 | 0.8686392 | NA | NA |
| Total | 31 | 12.030738 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.03491 0.034906 2.5203 999 0.136
Residuals 30 0.41549 0.013850
adonis2(neutral ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.78013 | 0.1607184 | 5.744856 | 0.001 |
| Residual | 30 | 9.29595 | 0.8392816 | NA | NA |
| Total | 31 | 11.07608 | 1.0000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.31336 0.313365 16.503 999 0.002 **
Residuals 30 0.56964 0.018988
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(phylo ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.572915 | 0.2890287 | 12.19579 | 0.001 |
| Residual | 30 | 3.869157 | 0.7109713 | NA | NA |
| Total | 31 | 5.442071 | 1.0000000 | NA | NA |
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.18493 0.184935 4.9256 999 0.037 *
Residuals 30 1.12638 0.037546
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(func ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(func))),
permutations = 999,
strata = subset_meta %>% arrange(match(Tube_code,labels(func))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.120331 | 0.3613445 | 16.97368 | 0.001 |
| Residual | 30 | 1.980120 | 0.6386555 | NA | NA |
| Total | 31 | 3.100451 | 1.0000000 | NA | NA |
4.3.3 Structural zeros
acclimation_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point == "1_Acclimation") %>%
dplyr::select(sample) %>%
pull()
antibiotics_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point == "2_Antibiotics") %>%
dplyr::select(sample) %>% pull()
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",subset_meta$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_acclimation = all(c_across(all_of(acclimation_samples)) == 0)) %>%
mutate(all_zeros_antibiotics = all(c_across(all_of(antibiotics_samples)) == 0)) %>%
mutate(average_acclimation = mean(c_across(all_of(acclimation_samples)), na.rm = TRUE)) %>%
mutate(average_antibiotics = mean(c_across(all_of(antibiotics_samples)), na.rm = TRUE)) %>%
filter(all_zeros_acclimation == TRUE || all_zeros_antibiotics==TRUE) %>%
mutate(present = case_when(
all_zeros_acclimation & !all_zeros_antibiotics ~ "antibiotics",
!all_zeros_acclimation & all_zeros_antibiotics ~ "acclimation",
!all_zeros_acclimation & !all_zeros_antibiotics ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "acclimation", average_acclimation, average_antibiotics)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Acclimation
structural_zeros %>%
filter(present=="acclimation") %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) # A tibble: 9 × 2
# Rowwise:
phylum sample_count
<chr> <int>
1 p__Bacillota_A 43
2 p__Bacteroidota 15
3 p__Bacillota 8
4 p__Pseudomonadota 7
5 p__Cyanobacteriota 3
6 p__Verrucomicrobiota 2
7 p__Bacillota_B 1
8 p__Bacillota_C 1
9 p__Fusobacteriota 1
Antibiotics
# A tibble: 4 × 7
# Rowwise:
genome phylum class order family genus species
<chr> <chr> <chr> <chr> <chr> <chr> <chr>
1 AH1_2nd_7:bin_000003 p__Pseudomonadota c__Gammaproteobacteria o__Enterobacterales f__Enterobacteriaceae g__Proteus s__Proteus cibarius
2 LI1_2nd_7:bin_000001 p__Bacillota_A c__Clostridia o__Clostridiales f__Clostridiaceae g__Sarcina s__
3 AH1_2nd_7:bin_000055 p__Bacillota c__Bacilli o__Mycoplasmatales f__Mycoplasmoidaceae g__Ureaplasma s__
4 AH1_2nd_13:bin_000025 p__Bacillota_A c__Clostridia o__Christensenellales f__UBA3700 g__ s__
Community elements differences in structural zeros
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)# A tibble: 0 × 3
# ℹ 3 variables: trait <chr>, p_value <dbl>, p_adjust <dbl>
4.3.4 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
filter(Population == "Cold_wet" & time_point %in% c("1_Acclimation", "2_Antibiotics") )%>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "time_point",
rand_formula = "(1|individual)",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output$res %>%
dplyr::select(taxon, lfc_time_point2_Antibiotics, p_time_point2_Antibiotics) %>%
filter(p_time_point2_Antibiotics < 0.05) %>%
dplyr::arrange(p_time_point2_Antibiotics) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_time_point2_Antibiotics)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_time_point2_Antibiotics, y=forcats::fct_reorder(genome,lfc_time_point2_Antibiotics), fill=phylum)) + #forcats::fct_rev()
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
# coord_flip()+
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))Phyla of the significant MAGs
ancombc_rand_table_mag%>%
filter(lfc_time_point2_Antibiotics<0) %>%
count(phylum, name = "sample_count") %>%
arrange(desc(sample_count)) phylum sample_count
1 Bacteroidota 14
2 Bacillota_A 4
3 Bacillota 2
4 Campylobacterota 1
phylum genus species
1 Campylobacterota Helicobacter_J
2 Bacillota_A
3 Bacteroidota Bacteroides
4 Bacteroidota Odoribacter
5 Bacteroidota Bacteroides
6 Bacillota Coprobacillus
7 Bacillota
8 Bacteroidota Phocaeicola
9 Bacteroidota Bacteroides
10 Bacteroidota Phocaeicola
11 Bacteroidota Odoribacter
12 Bacteroidota Bacteroides
13 Bacteroidota Bacteroides
14 Bacteroidota Parabacteroides
15 Bacteroidota
16 Bacteroidota Parabacteroides
17 Bacillota_A
18 Bacteroidota Alistipes
19 Bacillota_A
20 Bacillota_A
21 Bacteroidota Odoribacter
ancombc_rand_table_mag%>%
filter(lfc_time_point2_Antibiotics<0) %>%
count(genus, name = "sample_count") %>%
arrange(desc(sample_count)) genus sample_count
1 6
2 Bacteroides 5
3 Odoribacter 3
4 Parabacteroides 2
5 Phocaeicola 2
6 Alistipes 1
7 Coprobacillus 1
8 Helicobacter_J 1
4.3.5 Differences in the functional capacity
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
filter(Population == "Cold_wet" & treatment %in% c("Acclimation", "Antibiotics"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.3.5.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.331 0.0320
2 Antibiotics 0.244 0.134
MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94332, p-value = 0.09304
Wilcoxon rank sum exact test
data: value by treatment
W = 190, p-value = 0.01768
alternative hypothesis: true location shift is not equal to 0
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=Acclimation-Antibiotics)%>%
mutate(group_color = ifelse(Difference <0, "Antibiotics","Acclimation")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=treatment_colors) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")4.3.5.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 Acclimation 0.329 0.0320
2 Antibiotics 0.255 0.126
MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.95737, p-value = 0.2324
Wilcoxon rank sum exact test
data: value by treatment
W = 188, p-value = 0.02195
alternative hypothesis: true location shift is not equal to 0
function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))| Code_function | Acclimation | Antibiotics | Function | Difference | treatment |
|---|---|---|---|---|---|
| B06 | 0.68110280 | 0.47325490 | Organic anion biosynthesis | 0.20784790 | Acclimation |
| B02 | 0.59963040 | 0.41562100 | Amino acid biosynthesis | 0.18400940 | Acclimation |
| D02 | 0.38587770 | 0.20636760 | Polysaccharide degradation | 0.17951010 | Acclimation |
| S03 | 0.27103471 | 0.09357224 | Spore | 0.17746247 | Acclimation |
| B01 | 0.84215140 | 0.69078910 | Nucleic acid biosynthesis | 0.15136230 | Acclimation |
| B07 | 0.44558700 | 0.30432910 | Vitamin biosynthesis | 0.14125790 | Acclimation |
| B08 | 0.44596150 | 0.32044710 | Aromatic compound biosynthesis | 0.12551440 | Acclimation |
| D09 | 0.25533360 | 0.13438350 | Antibiotic degradation | 0.12095010 | Acclimation |
| B04 | 0.54457570 | 0.42921780 | SCFA biosynthesis | 0.11535790 | Acclimation |
| D03 | 0.29210750 | 0.20698550 | Sugar degradation | 0.08512200 | Acclimation |
| D05 | 0.28856080 | 0.22258820 | Amino acid degradation | 0.06597260 | Acclimation |
| B10 | 0.05947806 | 0.03621687 | Antibiotic biosynthesis | 0.02326119 | Acclimation |
4.4 Does antibiotic administration remove the differences found in the two populations?
4.4.1 Alpha diversity
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(time_point == "2_Antibiotics" )alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(time_point == "2_Antibiotics" ) %>%
ggplot(aes(y = value, x = Population, color=Population, fill=Population)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(values=treatment_colors)+
scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")4.4.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(time_point == "2_Antibiotics") %>%
select(Tube_code) %>%
pull()
subset_meta <- sample_metadata %>%
filter(time_point == "2_Antibiotics")Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$Population) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.015377 0.0153774 6.8934 999 0.011 *
Residuals 21 0.046845 0.0022307
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(richness ~ Population,
data = subset_meta %>% arrange(match(Tube_code,labels(richness))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.361743 | 0.1533845 | 3.80465 | 0.001 |
| Residual | 21 | 7.516224 | 0.8466155 | NA | NA |
| Total | 22 | 8.877966 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$Population) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.030649 0.0306488 3.8694 999 0.069 .
Residuals 21 0.166339 0.0079209
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(neutral ~ Population,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.787698 | 0.208772 | 5.541022 | 0.001 |
| Residual | 21 | 6.775221 | 0.791228 | NA | NA |
| Total | 22 | 8.562919 | 1.000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$Population) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.012165 0.012165 0.9979 999 0.354
Residuals 21 0.256012 0.012191
adonis2(phylo ~ Population,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.8970751 | 0.1890846 | 4.896659 | 0.001 |
| Residual | 21 | 3.8472307 | 0.8109154 | NA | NA |
| Total | 22 | 4.7443057 | 1.0000000 | NA | NA |
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$Population) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.01804 0.018036 0.4397 999 0.517
Residuals 21 0.86147 0.041023
adonis2(func ~ Population,
data = subset_meta %>% arrange(match(Tube_code,labels(func))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.0184103 | 0.009818749 | 0.2082384 | 0.678 |
| Residual | 21 | 1.8566048 | 0.990181251 | NA | NA |
| Total | 22 | 1.8750151 | 1.000000000 | NA | NA |
4.5 Are the microbial communities similar in both donor samples?
samples_to_keep <- sample_metadata %>%
filter(type =="Hot_control" & time_point %in% c( "3_Transplant1", "4_Transplant2"))%>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type =="Hot_control" & time_point %in% c( "3_Transplant1", "4_Transplant2"))alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(type =="Hot_control" & time_point %in% c( "3_Transplant1", "4_Transplant2")) %>%
ggplot(aes(y = value, x = time_point, color=time_point, fill=time_point)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(values=treatment_colors)+
scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00014 0.0001396 0.0173 999 0.912
Residuals 13 0.10481 0.0080621
adonis2(richness ~ time_point,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.06889268 | 0.03010712 | 0.4035421 | 0.981 |
| Residual | 13 | 2.21935895 | 0.96989288 | NA | NA |
| Total | 14 | 2.28825162 | 1.00000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000655 0.0006546 0.0514 999 0.816
Residuals 13 0.165409 0.0127237
adonis2(neutral ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(neutral))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.09294671 | 0.03882177 | 0.5250671 | 0.818 |
| Residual | 13 | 2.30124351 | 0.96117823 | NA | NA |
| Total | 14 | 2.39419022 | 1.00000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.003691 0.0036912 0.3071 999 0.673
Residuals 13 0.156266 0.0120205
adonis2(phylo ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(phylo))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.009773631 | 0.01921019 | 0.2546239 | 0.841 |
| Residual | 13 | 0.498999565 | 0.98078981 | NA | NA |
| Total | 14 | 0.508773196 | 1.00000000 | NA | NA |
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$time_point) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.013592 0.0135920 1.5644 999 0.195
Residuals 13 0.112945 0.0086881
adonis2(func ~ time_point,
data = subset_meta %>% arrange(match(Tube_code,labels(func))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | -0.003158048 | -0.0167584 | -0.2142684 | 0.821 |
| Residual | 13 | 0.191603701 | 1.0167584 | NA | NA |
| Total | 14 | 0.188445652 | 1.0000000 | NA | NA |
4.6 Does the donor sample maintain the microbial community found in acclimation?
4.6.1 Alpha diversity
sample_metadata <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Transplant1", "Transplant2") ~ "Transplant",
TRUE ~ treatment
))
samples_to_keep <- sample_metadata %>%
filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant"))alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(type == "Hot_control" & treatment %in% c("Acclimation","Transplant")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(values=treatment_colors)+
scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")4.6.2 Beta diversity
Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00836 0.0083600 1.4409 999 0.255
Residuals 22 0.12764 0.0058019
adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.2657351 | 0.06396199 | 1.503319 | 0.096 |
| Residual | 22 | 3.8888430 | 0.93603801 | NA | NA |
| Total | 23 | 4.1545781 | 1.00000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000661 0.0006614 0.0736 999 0.793
Residuals 22 0.197804 0.0089911
adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3164816 | 0.07596593 | 1.808646 | 0.07 |
| Residual | 22 | 3.8496178 | 0.92403407 | NA | NA |
| Total | 23 | 4.1660995 | 1.00000000 | NA | NA |
neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
# scale_color_manual(values = treatment_colors,labels=c("Cold_wet" = "Cold wet", "Hot_dry" = "Hot dry")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00204 0.0020405 0.2622 999 0.71
Residuals 22 0.17118 0.0077807
adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.0412056 | 0.056244 | 1.31111 | 0.268 |
| Residual | 22 | 0.6914166 | 0.943756 | NA | NA |
| Total | 23 | 0.7326222 | 1.000000 | NA | NA |
phylo %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000463 0.0004627 0.0465 999 0.839
Residuals 22 0.219125 0.0099602
adonis2(func ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(func))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.08931756 | 0.2172817 | 6.107175 | 0.056 |
| Residual | 22 | 0.32175047 | 0.7827183 | NA | NA |
| Total | 23 | 0.41106803 | 1.0000000 | NA | NA |
func %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
scale_color_manual(values = treatment_colors,labels=c("Cold_wet" = "Cold wet", "Hot_dry" = "Hot dry")) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)4.7 Is the donor sample microbiota different to recipient microbiota?
4.7.1 Alpha diversity
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant"))alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Transplant")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(values=treatment_colors)+
scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
stat_compare_means(method = "wilcox.test", show.legend = F, size = 3, label.x = c(1.5))+
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
# Increase plot size
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")4.7.2 Beta diversity
Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.12795 0.127946 13.179 999 0.001 ***
Residuals 20 0.19416 0.009708
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.777815 | 0.2760196 | 7.625059 | 0.001 |
| Residual | 20 | 4.663086 | 0.7239804 | NA | NA |
| Total | 21 | 6.440902 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.071502 0.071502 5.4967 999 0.031 *
Residuals 20 0.260166 0.013008
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 2.112572 | 0.3210813 | 9.458609 | 0.001 |
| Residual | 20 | 4.466983 | 0.6789187 | NA | NA |
| Total | 21 | 6.579556 | 1.0000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04731 0.047307 2.6157 999 0.115
Residuals 20 0.36172 0.018086
adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3778248 | 0.239656 | 6.303884 | 0.001 |
| Residual | 20 | 1.1987049 | 0.760344 | NA | NA |
| Total | 21 | 1.5765298 | 1.000000 | NA | NA |
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.000167 0.0001673 0.0139 999 0.901
Residuals 20 0.241060 0.0120530
adonis2(func ~ treatment,
data = subset_meta %>% arrange(match(Tube_code,labels(func))),
permutations = 999) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.07837549 | 0.1797296 | 4.382204 | 0.088 |
| Residual | 20 | 0.35769892 | 0.8202704 | NA | NA |
| Total | 21 | 0.43607441 | 1.0000000 | NA | NA |
4.8 Does FMT change the microbial community over time?
4.8.1 Alpha diversity
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Treatment" & treatment %in% c("Post-FMT1","Post-FMT2") ) alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(type=="Treatment" & treatment %in% c("Post-FMT1", "Post-FMT2" )) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
scale_color_manual(values=treatment_colors)+
scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")Richness
Model_rich <- glmer.nb(richness ~ treatment+(1|individual), data = alpha_div_meta)
summary(Model_rich)Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: Negative Binomial(4.3876) ( log )
Formula: richness ~ treatment + (1 | individual)
Data: alpha_div_meta
AIC BIC logLik deviance df.resid
171.0 174.3 -81.5 163.0 13
Scaled residuals:
Min 1Q Median 3Q Max
-1.73495 -0.71265 -0.07086 0.40744 1.88240
Random effects:
Groups Name Variance Std.Dev.
individual (Intercept) 1.073e-11 3.275e-06
Number of obs: 17, groups: individual, 9
Fixed effects:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.9343 0.1759 22.369 <2e-16 ***
treatmentPost-FMT2 0.4052 0.2402 1.687 0.0917 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
trtmnP-FMT2 -0.732
optimizer (Nelder_Mead) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')
Neutral
Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_neutral)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
128.7946 131.6268 -60.3973
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 3.497472 11.25384
Fixed effects: neutral ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 24.65947 4.164756 8 5.920988 0.0004
treatmentPost-FMT2 14.87494 5.482531 7 2.713151 0.0301
Correlation:
(Intr)
treatmentPost-FMT2 -0.7
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.76462323 -0.55701822 -0.04763166 0.55267588 1.30333436
Number of Observations: 17
Number of Groups: 9
Phylogenetic
Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_phylo)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
56.17121 59.00341 -24.0856
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.6979711 0.8383981
Fixed effects: phylogenetic ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 4.163606 0.3820859 8 10.897043 0.0000
treatmentPost-FMT2 0.916226 0.4122640 7 2.222426 0.0617
Correlation:
(Intr)
treatmentPost-FMT2 -0.583
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.1362905 -0.5779277 -0.1466831 0.4812394 2.3357705
Number of Observations: 17
Number of Groups: 9
Functional
Model_funct <- lme(fixed = functional ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
summary(Model_funct)Linear mixed-effects model fit by REML
Data: alpha_div_meta
AIC BIC logLik
-32.72222 -29.89002 20.36111
Random effects:
Formula: ~1 | individual
(Intercept) Residual
StdDev: 0.06160077 0.0305115
Fixed effects: functional ~ treatment
Value Std.Error DF t-value p-value
(Intercept) 1.4877804 0.02341751 8 63.53281 0.0000
treatmentPost-FMT2 0.0320561 0.01517203 7 2.11284 0.0725
Correlation:
(Intr)
treatmentPost-FMT2 -0.357
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-1.07361856 -0.66739366 -0.04576292 0.53118788 1.61981397
Number of Observations: 17
Number of Groups: 9
4.8.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Post-FMT1", "Post-FMT2")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Post-FMT1", "Post-FMT2"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.017610 0.017610 2.8225 999 0.102
Residuals 15 0.093585 0.006239
adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3560084 | 0.07968734 | 1.298809 | 0.00390625 |
| Residual | 15 | 4.1115570 | 0.92031266 | NA | NA |
| Total | 16 | 4.4675655 | 1.00000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.009828 0.0098278 0.7113 999 0.436
Residuals 15 0.207263 0.0138175
adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.3149954 | 0.08110541 | 1.323962 | 0.00390625 |
| Residual | 15 | 3.5687823 | 0.91889459 | NA | NA |
| Total | 16 | 3.8837776 | 1.00000000 | NA | NA |
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.010955 0.010955 0.6602 999 0.545
Residuals 15 0.248892 0.016593
adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.06391536 | 0.09449338 | 1.565312 | 0.0703125 |
| Residual | 15 | 0.61248501 | 0.90550662 | NA | NA |
| Total | 16 | 0.67640037 | 1.00000000 | NA | NA |
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00151 0.0015101 0.1075 999 0.727
Residuals 15 0.21063 0.0140417
adonis2(func ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(func))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(func))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.02861781 | 0.07980286 | 1.300855 | 0.3203125 |
| Residual | 15 | 0.32998847 | 0.92019714 | NA | NA |
| Total | 16 | 0.35860627 | 1.00000000 | NA | NA |
4.9 Do FMT change the microbial community compared to antibiotics and acclimation?
4.9.1 Alpha diversity
alpha_div_meta <- alpha_div %>%
left_join(sample_metadata, by = join_by(sample == sample))%>%
filter(type == "Treatment" & treatment %in% c("Antibiotics", "Acclimation","Post-FMT1") ) alpha_div %>%
pivot_longer(-sample, names_to = "metric", values_to = "value") %>%
left_join(., sample_metadata, by = "sample") %>%
mutate(metric=factor(metric,levels=c("richness","neutral","phylogenetic", "functional"))) %>%
filter(type=="Treatment" & treatment %in% c( "Antibiotics","Acclimation", "Post-FMT1")) %>%
ggplot(aes(y = value, x = treatment, color=treatment, fill=treatment)) +
geom_jitter(width = 0.2, show.legend = FALSE) +
geom_boxplot(width = 0.5, alpha=0.5, outlier.shape = NA, show.legend = FALSE) +
# scale_color_manual(values=treatment_colors)+
# scale_fill_manual(values=treatment_colors) +
facet_wrap(. ~ metric, scales = "free") +
theme_classic() +
theme(
strip.background = element_blank(),
panel.grid.minor.x = element_line(size = .1, color = "grey"),
axis.title.x = element_blank(),
axis.title.y = element_text(size=10),
axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(size = 10),
axis.text = element_text(size = 8),
axis.title = element_text(size = 8)
) +
ylab("Alpha diversity")Richness: Problems to calculate
Model_rich <- glmer.nb(richness ~ treatment + (1|individual), data = alpha_div_meta)
#summary(Model_rich)
emmeans(Model_rich, pairwise ~ treatment)Neutral
Model_neutral <- lme(fixed = neutral ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
#summary(Model_neutral)
emmeans(Model_neutral, pairwise ~ treatment)$emmeans
treatment emmean SE df lower.CL upper.CL
Acclimation 17.31 3.42 8 9.431 25.2
Antibiotics 8.44 3.86 8 -0.451 17.3
Post-FMT1 24.60 3.62 8 16.256 32.9
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Antibiotics 8.87 4.74 13 1.869 0.1868
Acclimation - (Post-FMT1) -7.29 4.55 13 -1.601 0.2800
Antibiotics - (Post-FMT1) -16.16 4.85 13 -3.333 0.0139
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates
Phylogenetic
Model_phylo <- lme(fixed = phylogenetic ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
#summary(Model_phylo)
emmeans(Model_phylo, pairwise ~ treatment)$emmeans
treatment emmean SE df lower.CL upper.CL
Acclimation 5.50 0.505 8 4.34 6.67
Antibiotics 4.30 0.571 8 2.98 5.62
Post-FMT1 4.15 0.535 8 2.91 5.38
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Antibiotics 1.203 0.716 13 1.680 0.2496
Acclimation - (Post-FMT1) 1.357 0.688 13 1.973 0.1583
Antibiotics - (Post-FMT1) 0.154 0.733 13 0.210 0.9760
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates
Functional
Model_funct <- lme(fixed = functional ~ treatment, data = alpha_div_meta,
random = ~ 1 | individual)
#summary(Model_funct)
emmeans(Model_funct, pairwise ~ treatment)$emmeans
treatment emmean SE df lower.CL upper.CL
Acclimation 1.49 0.0437 8 1.39 1.60
Antibiotics 1.47 0.0495 8 1.35 1.58
Post-FMT1 1.49 0.0463 8 1.38 1.60
Degrees-of-freedom method: containment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Acclimation - Antibiotics 0.02530 0.0660 13 0.383 0.9227
Acclimation - (Post-FMT1) 0.00618 0.0637 13 0.097 0.9948
Antibiotics - (Post-FMT1) -0.01913 0.0678 13 -0.282 0.9572
Degrees-of-freedom method: containment
P value adjustment: tukey method for comparing a family of 3 estimates
4.9.2 Beta diversity
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Antibiotics", "Post-FMT1")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("Acclimation", "Antibiotics", "Post-FMT1"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.01582 0.0079101 1.05 999 0.383
Residuals 21 0.15821 0.0075336
adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.840841 | 0.2018508 | 2.655435 | 0.001 |
| Residual | 21 | 7.278967 | 0.7981492 | NA | NA |
| Total | 23 | 9.119808 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.8287871 2.293722 0.1407734 0.002 0.006 *
2 Acclimation vs Post-FMT1 1 0.9572484 2.939042 0.1638350 0.001 0.003 *
3 Antibiotics vs Post-FMT1 1 0.9764237 2.751191 0.1746656 0.001 0.003 *
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.01796 0.0089778 0.5959 999 0.592
Residuals 21 0.31640 0.0150667
adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 2.357142 | 0.2699241 | 3.882065 | 0.001 |
| Residual | 21 | 6.375470 | 0.7300759 | NA | NA |
| Total | 23 | 8.732613 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.8814932 2.747044 0.1640316 0.001 0.003 *
2 Acclimation vs Post-FMT1 1 1.3133104 4.634354 0.2360329 0.001 0.003 *
3 Antibiotics vs Post-FMT1 1 1.3427497 4.355528 0.2509591 0.002 0.006 *
Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.15410 0.077052 2.6434 999 0.08 .
Residuals 21 0.61214 0.029149
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.327076 | 0.3647564 | 6.029092 | 0.001 |
| Residual | 21 | 2.311178 | 0.6352436 | NA | NA |
| Total | 23 | 3.638254 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 0.6885926 5.124832 0.2679674 0.005 0.015 .
2 Acclimation vs Post-FMT1 1 0.2548027 3.292748 0.1800029 0.045 0.135
3 Antibiotics vs Post-FMT1 1 1.1000473 9.048069 0.4103792 0.002 0.006 *
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 2 0.23408 0.117038 4.3947 999 0.024 *
Residuals 21 0.55927 0.026632
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(func ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(func))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(func))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 2 | 1.636530 | 0.5730275 | 14.09175 | 0.001 |
| Residual | 21 | 1.219406 | 0.4269725 | NA | NA |
| Total | 23 | 2.855936 | 1.0000000 | NA | NA |
pairs Df SumsOfSqs F.Model R2 p.value p.adjusted sig
1 Acclimation vs Antibiotics 1 1.11236925 14.811924 0.5140901 0.003 0.009 *
2 Acclimation vs Post-FMT1 1 0.05914673 2.630593 0.1492061 0.166 0.498
3 Antibiotics vs Post-FMT1 1 1.36488774 16.896102 0.5651607 0.001 0.003 *
4.10 Is the gut microbiota similar to the donor after FMT?
4.10.1 Beta diversity
sample_metadata <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
TRUE ~ treatment
))
samples_to_keep <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant")) %>%
select(sample) %>%
pull()
subset_meta <- sample_metadata %>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))Richness
richness <- as.matrix(beta_q0n$S)
richness <- as.dist(richness[rownames(richness) %in% samples_to_keep,
colnames(richness) %in% samples_to_keep])
betadisper(richness, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.11100 0.111002 11.752 999 0.001 ***
Residuals 28 0.26447 0.009445
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(richness ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(richness))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(richness))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.184578 | 0.1548451 | 5.13002 | 0.001 |
| Residual | 28 | 6.465508 | 0.8451549 | NA | NA |
| Total | 29 | 7.650086 | 1.0000000 | NA | NA |
Neutral
neutral <- as.matrix(beta_q1n$S)
neutral <- as.dist(neutral[rownames(neutral) %in% samples_to_keep,
colnames(neutral) %in% samples_to_keep])
betadisper(neutral, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.04408 0.044082 3.1744 999 0.092 .
Residuals 28 0.38883 0.013887
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis2(neutral ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(neutral))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(neutral))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 1.147077 | 0.1608783 | 5.368223 | 0.001 |
| Residual | 28 | 5.983012 | 0.8391217 | NA | NA |
| Total | 29 | 7.130089 | 1.0000000 | NA | NA |
neutral %>%
vegan::metaMDS(., trymax = 500, k = 2, trace=0) %>%
vegan::scores() %>%
as_tibble(., rownames = "sample") %>%
dplyr::left_join(subset_meta, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
mutate(x_cen = mean(NMDS1, na.rm = TRUE)) %>%
mutate(y_cen = mean(NMDS2, na.rm = TRUE)) %>%
ungroup() %>%
ggplot(aes(x = NMDS1, y = NMDS2, color = treatment)) +
geom_point(size = 4) +
# scale_color_manual(values = treatment_colors) +
geom_segment(aes(x = x_cen, y = y_cen, xend = NMDS1, yend = NMDS2), alpha = 0.9, show.legend = FALSE) +
theme(
axis.title = element_text(size = 12, face = "bold"),
axis.text = element_text(size = 10),
panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
legend.text = element_text(size = 10),
legend.title = element_text(size = 12),
legend.position = "right", legend.box = "vertical"
)Phylogenetic
phylo <- as.matrix(beta_q1p$S)
phylo <- as.dist(phylo[rownames(phylo) %in% samples_to_keep,
colnames(phylo) %in% samples_to_keep])
betadisper(phylo, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00007 0.0000689 0.0044 999 0.948
Residuals 28 0.43637 0.0155847
adonis2(phylo ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(phylo))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(phylo))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.1298187 | 0.1018776 | 3.17615 | 0.001 |
| Residual | 28 | 1.1444432 | 0.8981224 | NA | NA |
| Total | 29 | 1.2742619 | 1.0000000 | NA | NA |
Functional
func <- as.matrix(beta_q1f$S)
func <- as.dist(func[rownames(func) %in% samples_to_keep,
colnames(func) %in% samples_to_keep])
betadisper(func, subset_meta$treatment) %>% permutest(.)
Permutation test for homogeneity of multivariate dispersions
Permutation: free
Number of permutations: 999
Response: Distances
Df Sum Sq Mean Sq F N.Perm Pr(>F)
Groups 1 0.00121 0.0012109 0.1039 999 0.78
Residuals 28 0.32638 0.0116565
adonis2(func ~ treatment,
data = subset_meta %>% arrange(match(sample,labels(func))),
permutations = 999,
strata = subset_meta %>% arrange(match(sample,labels(func))) %>% pull(individual)) %>%
broom::tidy() %>%
tt()| term | df | SumOfSqs | R2 | statistic | p.value |
|---|---|---|---|---|---|
| Model | 1 | 0.01515762 | 0.02696066 | 0.7758149 | 0.481 |
| Residual | 28 | 0.54705472 | 0.97303934 | NA | NA |
| Total | 29 | 0.56221234 | 1.00000000 | NA | NA |
4.10.2 Structural zeros
FMT_samples <- sample_metadata %>%
filter(type == "Treatment" & treatment == "FMT") %>%
dplyr::select(sample) %>%
pull()
Transplant_samples <- sample_metadata %>%
filter(type == "Treatment" & treatment =="Transplant") %>%
dplyr::select(sample) %>% pull()
structural_zeros <- genome_counts_filt %>%
select(one_of(c("genome",subset_meta$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
rowwise() %>% #compute for each row (genome)
mutate(all_zeros_FMT = all(c_across(all_of(FMT_samples)) == 0)) %>%
mutate(all_zeros_Transplant = all(c_across(all_of(Transplant_samples)) == 0)) %>%
mutate(average_FMT = mean(c_across(all_of(FMT_samples)), na.rm = TRUE)) %>%
mutate(average_Transplant = mean(c_across(all_of(Transplant_samples)), na.rm = TRUE)) %>%
filter(all_zeros_FMT == TRUE || all_zeros_Transplant==TRUE) %>%
mutate(present = case_when(
all_zeros_FMT & !all_zeros_Transplant ~ "Transplant",
!all_zeros_FMT & all_zeros_Transplant ~ "FMT",
!all_zeros_FMT & !all_zeros_Transplant ~ "None",
TRUE ~ NA_character_
)) %>%
mutate(average = ifelse(present == "FMT", average_FMT, average_Transplant)) %>%
dplyr::select(genome, present, average) %>%
left_join(genome_metadata, by=join_by(genome==genome)) %>%
arrange(present,-average)Structural zeros in each group
fmt <- structural_zeros %>%
filter(present=="FMT") %>%
count(phylum, name = "FMT") %>%
arrange(desc(FMT))
fmt_f <- structural_zeros %>%
filter(present=="FMT") %>%
count(family, name = "FMT") %>%
arrange(desc(FMT)) structural_zeros %>%
filter(present=="Transplant") %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt, by="phylum" ) %>%
mutate(across(everything(), ~ replace_na(., 0))) %>%
as.data.frame() %>%
summarise(across(where(is.numeric), ~ sum(.x, na.rm = TRUE))) Transplant FMT
1 52 55
Phyla to which the structural zeros belong in each group
structural_zeros %>%
filter(present=="Transplant") %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt, by="phylum" ) %>%
mutate(across(everything(), ~ replace_na(., 0))) %>%
tt()| phylum | Transplant | FMT |
|---|---|---|
| p__Bacillota_A | 20 | 21 |
| p__Bacillota | 16 | 10 |
| p__Pseudomonadota | 6 | 11 |
| p__Bacteroidota | 3 | 6 |
| p__Desulfobacterota | 2 | 2 |
| p__Bacillota_B | 1 | 0 |
| p__Chlamydiota | 1 | 0 |
| p__Cyanobacteriota | 1 | 0 |
| p__Fusobacteriota | 1 | 0 |
| p__Verrucomicrobiota | 1 | 2 |
| p__Actinomycetota | 0 | 1 |
| p__Bacillota_C | 0 | 1 |
| p__Campylobacterota | 0 | 1 |
Families to which the structural zeros belong in each group
structural_zeros %>%
filter(present=="Transplant") %>%
count(family, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt_f, by="family" ) %>%
mutate(across(everything(), ~ replace_na(., 0))) %>%
tt()| family | Transplant | FMT |
|---|---|---|
| f__Lachnospiraceae | 7 | 9 |
| f__Erysipelotrichaceae | 6 | 5 |
| f__UBA660 | 6 | 0 |
| f__Enterobacteriaceae | 5 | 2 |
| f__CAG-508 | 3 | 0 |
| f__Ruminococcaceae | 3 | 3 |
| f__Anaerovoracaceae | 2 | 0 |
| f__Coprobacillaceae | 2 | 0 |
| f__Desulfovibrionaceae | 2 | 2 |
| f__Oscillospiraceae | 2 | 1 |
| f__Tannerellaceae | 2 | 1 |
| f__UBA1242 | 2 | 0 |
| f__ | 1 | 3 |
| f__Akkermansiaceae | 1 | 0 |
| f__CAG-239 | 1 | 2 |
| f__Chlamydiaceae | 1 | 0 |
| f__Enterococcaceae | 1 | 3 |
| f__Fusobacteriaceae | 1 | 0 |
| f__Gastranaerophilaceae | 1 | 0 |
| f__Marinifilaceae | 1 | 0 |
| f__Mycoplasmoidaceae | 1 | 1 |
| f__Peptococcaceae | 1 | 0 |
| f__Anaerotignaceae | 0 | 2 |
| f__Bacteroidaceae | 0 | 2 |
| f__LL51 | 0 | 2 |
| f__UBA3700 | 0 | 2 |
| f__Acutalibacteraceae | 0 | 1 |
| f__Arcobacteraceae | 0 | 1 |
| f__CAG-274 | 0 | 1 |
| f__CALVMC01 | 0 | 1 |
| f__Devosiaceae | 0 | 1 |
| f__Mycobacteriaceae | 0 | 1 |
| f__Pumilibacteraceae | 0 | 1 |
| f__RUG11792 | 0 | 1 |
| f__Rhizobiaceae | 0 | 1 |
| f__Rikenellaceae | 0 | 1 |
| f__Sphingobacteriaceae | 0 | 1 |
| f__Streptococcaceae | 0 | 1 |
| f__UBA1997 | 0 | 1 |
| f__UBA3830 | 0 | 1 |
| f__Weeksellaceae | 0 | 1 |
4.10.2.1 Functional capacities of the structural zeros
#Aggregate bundle-level GIFTs into the compound level
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% structural_zeros$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
#Aggregate element-level GIFTs into the function level
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
#Aggregate function-level GIFTs into overall Biosynthesis, Degradation and Structural GIFTs
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
TRUE ~ treatment
))%>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
filter(genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=FMT-Transplant)%>%
mutate(group_color = ifelse(Difference <0, "Transplant","FMT")) difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=treatment_colors) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Elevation")4.10.3 Differential abundance analysis: composition
phylo_samples <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
TRUE ~ treatment
))%>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant")) %>%
column_to_rownames("sample") %>%
sample_data()
phylo_genome <- genome_counts_filt %>%
filter(!genome %in% structural_zeros$genome) %>%
select(one_of(c("genome",rownames(phylo_samples)))) %>%
filter(rowSums(across(where(is.numeric), ~ . != 0)) > 0) %>%
select(genome, where(~ is.numeric(.) && sum(.) > 0)) %>%
column_to_rownames("genome") %>%
otu_table(., taxa_are_rows = TRUE)
phylo_taxonomy <- genome_metadata %>%
filter(genome %in% rownames(phylo_genome)) %>%
column_to_rownames("genome") %>%
dplyr::select(domain,phylum,class,order,family,genus,species) %>%
as.matrix() %>%
tax_table()
physeq_genome_filtered <- phyloseq(phylo_genome, phylo_taxonomy, phylo_samples)ancom_rand_output = ancombc2(data = physeq_genome_filtered,
assay_name = "counts",
tax_level = NULL,
fix_formula = "treatment",
p_adj_method = "holm",
pseudo_sens = TRUE,
prv_cut =0,
lib_cut = 0,
s0_perc = 0.05,
group = NULL,
struc_zero = FALSE,
neg_lb = FALSE,
alpha = 0.05,
n_cl = 2,
verbose = TRUE,
global = FALSE,
pairwise = FALSE,
dunnet = FALSE,
trend = FALSE,
iter_control = list(tol = 1e-5, max_iter = 20, verbose = FALSE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
mdfdr_control = list(fwer_ctrl_method = "holm", B = 100),
trend_control = NULL)taxonomy <- data.frame(physeq_genome_filtered@tax_table) %>%
rownames_to_column(., "taxon") %>%
mutate_at(vars(order, phylum, family, genus, species), ~ str_replace(., "[dpcofgs]__", ""))
ancombc_rand_table_mag <- ancom_rand_output$res %>%
dplyr::select(taxon, lfc_treatmentTransplant, p_treatmentTransplant) %>%
filter(p_treatmentTransplant < 0.05) %>%
dplyr::arrange(p_treatmentTransplant) %>%
merge(., taxonomy, by="taxon") %>%
mutate_at(vars(phylum, species), ~ str_replace(., "[dpcofgs]__", ""))%>%
dplyr::arrange(lfc_treatmentTransplant)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate_at(vars(phylum), ~ str_replace(., "[dpcofgs]__", "")) %>%
right_join(taxonomy, by=join_by(phylum == phylum)) %>%
dplyr::select(phylum, colors) %>%
mutate(colors = str_c(colors, "80")) %>% #add 80% alpha
unique() %>%
dplyr::arrange(phylum)
tax_table <- as.data.frame(unique(ancombc_rand_table_mag$phylum))
colnames(tax_table)[1] <- "phylum"
tax_color <- merge(tax_table, colors_alphabetic, by="phylum")%>%
dplyr::arrange(phylum) %>%
dplyr::select(colors) %>%
pull()Differential abundance MAGs in each group
fmt_count <- ancombc_rand_table_mag %>%
filter(lfc_treatmentTransplant<0) %>%
count(phylum, name = "FMT") %>%
arrange(desc(FMT))
ancombc_rand_table_mag %>%
filter(lfc_treatmentTransplant>0) %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt_count, by="phylum")%>%
mutate(across(where(is.numeric), ~ replace_na(., 0))) phylum Transplant FMT
1 Bacteroidota 13 5
2 Bacillota_A 4 17
3 Bacillota 3 1
4 Campylobacterota 1 1
5 Desulfobacterota 1 2
6 Spirochaetota 1 0
7 Pseudomonadota 0 5
8 Cyanobacteriota 0 2
9 Bacillota_B 0 1
10 Bacillota_C 0 1
ancombc_rand_table_mag %>%
filter(lfc_treatmentTransplant>0) %>%
count(phylum, name = "Transplant") %>%
arrange(desc(Transplant)) %>%
full_join(.,fmt_count, by="phylum")%>%
mutate(across(where(is.numeric), ~ replace_na(., 0))) %>%
as.data.frame() %>%
summarise(across(where(is.numeric), ~ sum(.x, na.rm = TRUE))) Transplant FMT
1 23 35
ancombc_rand_table_mag%>%
mutate(genome=factor(taxon,levels=ancombc_rand_table_mag$taxon)) %>%
ggplot(., aes(x=lfc_treatmentTransplant, y=forcats::fct_reorder(genome,lfc_treatmentTransplant), fill=phylum)) + #forcats::fct_rev()
geom_col() +
scale_fill_manual(values=tax_color) +
geom_hline(yintercept=0) +
# coord_flip()+
theme(panel.background = element_blank(),
axis.line = element_line(size = 0.5, linetype = "solid", colour = "black"),
axis.text.x = element_text(size = 12),
axis.text.y = element_text(size = 8),
axis.title = element_text(size = 14, face = "bold"),
legend.text = element_text(size = 12),
legend.title = element_text(size = 14, face = "bold"),
legend.position = "right", legend.box = "vertical")+
xlab("log2FoldChange") +
ylab("Species")+
guides(fill=guide_legend(title="Phylum"))4.10.4 Differences in the functional capacities
GIFTs_elements <- to.elements(genome_gifts,GIFT_db)
GIFTs_elements_filtered <- GIFTs_elements[rownames(GIFTs_elements) %in% genome_counts_filt$genome,]
GIFTs_elements_filtered <- as.data.frame(GIFTs_elements_filtered) %>%
select_if(~ !is.numeric(.) || sum(.) != 0)
GIFTs_functions <- to.functions(GIFTs_elements_filtered,GIFT_db)
GIFTs_domains <- to.domains(GIFTs_functions,GIFT_db)
sample_sub <- sample_metadata %>%
mutate(treatment=case_when(
treatment %in% c("Post-FMT1", "Post-FMT2") ~ "FMT",
TRUE ~ treatment
))%>%
filter(type == "Treatment" & treatment %in% c("FMT", "Transplant"))
genome_counts_row <- genome_counts_filt %>%
mutate_at(vars(-genome),~./sum(.)) %>%
select(one_of(c("genome",sample_sub$sample))) %>%
column_to_rownames(., "genome")
GIFTs_elements_community <- to.community(GIFTs_elements_filtered,genome_counts_row,GIFT_db)
GIFTs_functions_community <- to.community(GIFTs_functions,genome_counts_row,GIFT_db)
GIFTs_domains_community <- to.community(GIFTs_domains,genome_counts_row,GIFT_db)4.10.4.1 Community elements differences:
MCI
GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 FMT 0.354 0.0255
2 Transplant 0.351 0.0426
MCI <- GIFTs_elements_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.94871, p-value = 0.1561
Wilcoxon rank sum exact test
data: value by treatment
W = 128, p-value = 0.4826
alternative hypothesis: true location shift is not equal to 0
element_gift <- GIFTs_elements_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
left_join(., sample_metadata[c(12,13)], by=join_by("sample"=="sample"))uniqueGIFT_db<- unique(GIFT_db[c(2,4,5,6)]) %>% unite("Function",Function:Element, sep= "_", remove=FALSE)
significant_elements <- element_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3)],by = join_by(trait == Code_element))
element_gift_t <- element_gift %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "trait")
element_gift_filt <- subset(element_gift_t, trait %in% significant_elements$trait) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "sample")%>%
left_join(., sample_metadata[c(12,13)], by = join_by(sample == sample))
difference_table <- element_gift_filt %>%
select(-sample) %>%
group_by(treatment) %>%
summarise(across(everything(), mean)) %>%
t() %>%
row_to_names(row_number = 1) %>%
as.data.frame() %>%
mutate_if(is.character, as.numeric) %>%
rownames_to_column(., "Elements") %>%
left_join(.,uniqueGIFT_db[c(1,3,4)],by = join_by(Elements == Code_element)) %>%
arrange(Function) %>%
mutate(Difference=FMT-Transplant)%>%
mutate(group_color = ifelse(Difference <0, "Transplant","FMT")) means_gift <- element_gift_filt %>%
select(-treatment) %>%
pivot_longer(!sample, names_to = "Elements", values_to = "abundance") %>%
left_join(sample_metadata, by=join_by(sample==sample)) %>%
group_by(treatment, Elements) %>%
summarise(mean=mean(abundance))
log_fold <- means_gift %>%
group_by(Elements) %>%
summarise(
logfc_fmt_transplant = log2(mean[treatment == "FMT"] / mean[treatment == "Transplant"])
) %>%
left_join(., difference_table, by="Elements")difference_table %>%
ggplot(aes(x=forcats::fct_reorder(Function,Difference), y=Difference, fill=group_color)) +
geom_col() +
# geom_point(size=4) +
scale_fill_manual(values=treatment_colors) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Mean difference")+
labs(fill="Treatment")log_fold %>%
filter(Elements!="D0611") %>%
ggplot(aes(x=forcats::fct_reorder(Function,logfc_fmt_transplant), y=logfc_fmt_transplant, fill=group_color)) +
geom_col() +
scale_fill_manual(values=treatment_colors) +
geom_hline(yintercept=0) +
coord_flip()+
theme(axis.text = element_text(size = 10),
axis.title = element_text(size = 12),
legend.position = "right",
panel.background = element_blank(),
panel.grid.major = element_line(size = 0.15, linetype = 'solid',
colour = "grey"))+
xlab("Microbial Functional Capacity") +
ylab("Log_fold")+
labs(fill="Treatment")4.10.4.2 Community functions differences
MCI
GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample)) %>%
group_by(treatment) %>%
summarise(MCI = mean(value), sd = sd(value))# A tibble: 2 × 3
treatment MCI sd
<chr> <dbl> <dbl>
1 FMT 0.349 0.0221
2 Transplant 0.354 0.0375
MCI <- GIFTs_functions_community %>%
rowMeans() %>%
as_tibble(., rownames = "sample") %>%
left_join(sample_metadata, by = join_by(sample == sample))
shapiro.test(MCI$value)
Shapiro-Wilk normality test
data: MCI$value
W = 0.87207, p-value = 0.001864
Wilcoxon rank sum exact test
data: value by treatment
W = 120, p-value = 0.7106
alternative hypothesis: true location shift is not equal to 0
function_gift <- GIFTs_functions_community %>%
as.data.frame() %>%
rownames_to_column(., "sample") %>%
merge(., sample_metadata[c(12,13)], by="sample")unique_funct_db<- GIFT_db[c(3,4,5)] %>%
distinct(Code_function, .keep_all = TRUE)
significant_functional <- function_gift %>%
pivot_longer(-c(sample,treatment), names_to = "trait", values_to = "value") %>%
group_by(trait) %>%
summarise(p_value = wilcox.test(value ~ treatment)$p.value) %>%
mutate(p_adjust=p.adjust(p_value, method="BH")) %>%
filter(p_adjust < 0.05)%>%
left_join(.,unique_funct_db[c(1,3)],by = join_by(trait == Code_function))